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The  time  dependent  mixed-layer  depth  and  temperature  re- 
sponse of  a  two-layer  hydrothermodynamic  ocean  model  to  a 
moving  hurricane  model  was  investigated.   The  hurricane  model 
was  that  of  Elsberry,  Pearson,  and  Corgnati  (1974).   The 
maximum  eye  radius  and  hurricane  intensity  was  bounded  to 
isolate  ocean  reaction  from  hurricane  variance.   The  ocean 
model  was  based  on  the  upwelling  model  of  O'Brien  and  Hurl- 
burt  (1972)  in  two  dimensions  and  was  located  along  the  path 
of  the  storm.   The  momentum  and  depth  equations  in  this 
model  are  treated  semi-implicit ly .   Mechanical  and  convective 
mixing  was  included  and  was  found  to  be  the  dominate  cause 
of  deepening  of  the  mixed  layer  and  cooling  of  the  ocean  sur- 
face.  Alternating  regions  of  upwelling  and  downwelling  were 
produced  in  the  wake  of  the  storm  as  a  result  of  inertial 
waves  caused  by  the  passage  of  the  storm.   The  cases  studied 
indicated  a  linear  relation  between  the  wavelength  and  the 
storm  translation  speed. 
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I.   INTRODUCTION 

It  is  generally  accepted  that  sea-surface  temperatures 
may  be  strongly  affected  by  tropical  storm  activity,  and 
vice  versa.   Sea-surface  temperatures  in  excess  of  26°C  are 
apparently  required  to  produce  hurricane-force  winds.   It 
has  also  been  observed  that  cooling  of  as  much  as  6°C  may  be 
present  in  the  wake  of  a  hurricane. 

One  of  the  first  systematic  efforts  to  investigate  the 
cooling  effect  was  done  after  the  passage  of  hurricane  Hilda, 
1964.   From  this  data,  Leipper  (1967)  found  clear  evidence 
supporting  previous  scattered  observations  which  indicated 
lowering  of  sea-surface  temperatures  in  the  wake  of  a  hurri- 
cane.  The  sea- surf ace  temperature  was  observed  to  drop  5  to 
6°C  from  the  temperature  present  before  the  storm  passage. 
The  observations  indicated  upwelling  from  about  60  m  in  depth, 
Associated  with  this  was  an  outward  transport  of  water  which 
occurred  in  the  warm  surface  layers  with  cooling  and  mixing 
of  these  layers  as  they  moved  outward.   Finally  this  outward 
moving  water  experienced  convergence  and  downwelling  around 
the  hurricane  area.   The  observations  were  taken  along  paths 
perpendicular  to  the  path  followed  by  the  storm.   This  would 
not  allow  certain  features  along  the  path  of  the  storm  to  be 
observed,  for  example,  the  internal  waves  predicted  by 
Geisler  (1970). 
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In  conjunction  with  project  Stormfury,  Black  and  Mallin- 
ger  (1972)  reported  on  data  taken  before,  during,  and  after 
the  passage  of  hurricane  Ginger,  1971.   These  oceanic  mea- 
surements were  taken  from  a  Navy  WP-3  aircraft  using  air- 
borne expendable  bathythermographs,  (AXBTs).   The  primary 
information  gained,  of  interest  in  this  thesis,  was  that  as 
the  storm  moved  faster  less  surface  cooling  occurred.   In 
fact,  there  appeared  to  be  a  cut-off  speed  related  to  the 
ratio  of  the  baroclinic  wave  speed  and  the  speed  of  propaga- 
tion of  the  storm.   Geisler  (1970)  predicted  this  ratio  to 
be  /2  whereas  the  data  obtained  suggested  a  ratio  of  3.   At 
the  slower  speeds  strong  upwelling  near  the  center  should 
occur.   While,  at  the  faster  speeds  the  upwelling  intensity 
decreases,  and  therefore  does  not  contribute  greatly  to  the 
surface  cooling. 

Sheets  (1974)  presented  evidence  indicating  internal 
gravity  waves  at  the  thermocline  in  the  wake  of  a  rapidly 
moving  hurricane,  Ellen,  1973.   A  total  of  18  AXBTs  were 
dropped  at  intervals  of  five  to  fifteen  nautical  miles.   In 
addition  to  the  internal  wave  indication,  a  surface  cooling 
of  .  5°C  to  1.5°C  was  observed.   The  amount  of  surface  cool- 
ing is  in  agreement  with  the  concept  that  a  rapidly  moving 
storm  does  not  cool  the  ocean  surface  to  the  same  extent  as 
slower  moving  storms.   In  the  same  paper,  however,  data  from 
hurricane  Irah,  a  slow  moving  storm  (7-10  kt ) ,  showed  nomi- 
nal cooling,  less  than  .5°C,  instead  of  the  greater  antici- 
pated cooling.   These  data  illustrate  that  many  factors  play 
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a  part  in  the  response  of  the  ocean  surface  temperature  to 
the  passage  of  a  hurricane. 

O'Brien  and  Reid  (1967)  developed  a  model  of  the  time- 
dependent  response  of  a  two-layer  (one  active  layer)  ocean 
model  under  a  stationary,  axially-symmetric  atmospheric 
vortex.   The  dynamic  response  of  the  ocean  was  also  assumed 
to  be  axially-symmetric,  and  therefore  solutions  in  two  in- 
dependent variables,  radial  distance  and  time,  could  be  ob- 
tained.  Gradient  winds  were  obtained  from  a  pressure 
distribution  function  resembling  actual  data.   The  surface 
stress  was  derived  using  these  winds.   The  results  of  this 
model  indicated  strong  upwelling  near  the  center  due  to  ve- 
locity divergence.   This  upwelling  was  strongest  near  the 
region  of  maximum  surface  stress.   After  48  hours  the  depth 
of  the  upper  layer  had  been  reduced  to  about  10  meters  in 
the  region  of  maximum  upwelling. 

Geisler  (1970)  developed  a  model  of  the  linear  response 
of  a  two-layer  ocean  to  a  stress  pattern  representing  a 
hurricane  moving  at  a  constant  speed.   When  the  speed  of 
movement  was  greater  than  the  baroclinic  wave  speed,  this 
model  produced  internal  waves  at  the  interface  of  the  mixed 
layer  and  the  thermocline  in  the  wake  of  the  storm.   The 
wavelength  of  these  disturbances  increased  as  the  speed  of 
movement  of  the  storm  increased.   The  disturbance  propagated 
in  the  same  direction  and  at  the  same  speed  as  the  storm. 
Another  feature  of  this  model  was  that  a  baroclinic  ridge 
was  left  in  the  ocean  in  the  wake  of  the  storm.   Due  to 
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upwelling  near  the  center  and  downwelling  in  the  surrounding 
area,  a  wake  of  reduced  mixed-layer  depth  near  the  center  of 
the  path  was  produced.   The  amplitude  of  this  ridge  was 
greater  with  decreasing  speed  of  movement  of  the  hurricane. 

The  models  discussed  thus  far  treated  only  the  response 
of  the  ocean  to  a  given  wind  stress.   In  an  effort  to  devel- 
op a  model  to  describe  the  interaction  of  the  ocean  and  the 
atmosphere  during  a  hurricane,  Elsberry,  Pearson  and  Corgnati 
(1974)  developed  an  axisymmetric  hurricane  model.   As  the 
atmospheric  portion  of  this  model  is  used  in  this  thesis  it 
will  be  described  later.   The  ocean  used  in  this  case  was  at 
a  fixed  temperature. 

Fraim  (1973)  introduced  a  mixed-layer  model  capable  of 
producing  time-dependent  values  of  sea-surface  temperature 
and  mixed-layer  depth  in  response  to  atmospheric  forcing. 
It  was  based  on  the  mixed-layer  models  of  Kraus  and  Turner 
(1967),  and  Denman  (1973),  which  predicted  the  mixed-layer 
temperature  and  depth  in  response  to  wind  mixing  and  convec- 
tion.  The  radial  current  was  defined  by  the  Ekman  current 
response  to  the  stress  at  the  surface.   In  this  model  the 
stress  input  to  wind  mixing,  t  ,  was  a  function  of  the  total 
stress,  t  ,  the  roughness  length,  Z0,  the  mixed-layer  depth, 
ho,  and  an  adjustable  constant  C. 


T.  ■  Ta  EXP(-  SZ7> 


The  excess  stress,  t   -  t  ,  was  used  to  produce  currents  from 

'   a     o  '  * 

Ekman  balance  considerations.   This  yielded  a  depth-averaged 


14 


radial  current,  and  due  to  symmetry,  this  is  the  only  cur- 
rent which  could  produce  upwelling.   Application  of  the 
continuity  equation  to  these  values  of  current  gave  the  up- 
welling  at  the  bottom  of  the  mixed  layer.   This  model  was 
used  to  study  the  effects  of  the  three  cooling  processes  in- 
volved in  lowering  the  sea  surface  temperatures.   These  are, 
entrainment  mixing  across  a  stable  layer,  convective  mixing, 
and  vertical  and  horizontal  advection.   The  results  of  these 
experiments  showed  that  both  advection  and  mixing  have  sig- 
nificant effects  on  the  sea  surface  temperature  and  mixed 
layer  depths,  and  as  expected  these  effects  are  most  pro- 
nounced in  the  region  of  maximum  surface  stress. 

Trapnell  (1974)  introduced  more  refinements  into  the 
Fraim  model,  the  most  notable  of  which  was  the  addition  of  a 
return  flow  in  the  lower  layer  and  the  ability  to  let  the 
surface  stress  move  across  the  ocean.   The  flow  in  the  lower 
layer  was  assumed  to  have  the  same  mass  transport  as  the 
upper  layer.   This  added  a  new  method  by  which  the  below 
layer  gradient  could  be  changed.   In  addition  to  upwelling 
and  downwelling,  there  now  existed  horizontal  advection. 
Movement  of  the  stress  and  cooling  functions  was  simulated 
by  translating  the  center  of  the  atmosphere  grid  relative  to 
the  ocean  grid.   This  adjustment  was  made  at  certain  time 
increments  to  simulate  different  speeds.   The  results  from 
Trapnell' s  experiments  suggest  that  entrainment  mixing  is 
the  primary  mechanisms  for  cooling  the  oceans,  except  within 
the  region  of  maximum  winds  where  intense  upwelling  becomes 
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important.   Also  it  is  seen  that  these  effects  are  strongly 
associated  with  slower  moving  storms.   Due  to  the  lack  of 
dynamics  this  model  did  not  exhibit  the  internal  waves  pre- 
dicted by  Geisler. 

One  of  the  primary  constraints  in  both  Fraim's  and 
Trapnell's  models  was  that  the  currents  used  were  derived 
from  Ekman  balance  to  a  given  surface  stress.   As  a  further 
refinement,  this  thesis  describes  a  two-layer  hydrothermo- 
dynamic  model  and  its  response  to  surface  stress  and  cooling 
representing  a  moving  hurricane. 
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II.   MODEL 

A.   OCEAN  MODEL 

The  dynamical  equations  in  the  ocean  model  and  the  meth- 
od of  solution  described  in  this  thesis  follow  the  coastal 
upwelling  model  of  O'Brien  and  Hurlburt  (1972).   However  the 
present  model  contains  a  term  to  include  thermal  and  mechan- 
ical mixing  due  to  upward  heat  flux  and  mechanical  energy 
input  at  the  surface. 

S    I 

dUi   .     dUi.    9  / ,   .,   \     r>     ,T   — T     .  .  9  U  1  /  t  \ 

jf- +  Ui9ir  +  ^(h^+h2>  =  fv>  +     phl     +A¥x^  (1) 


I    B 
x   x 


9u?  .    9u?    3  ,,     , ,     N 

J?  +  U2lt   +g^(hi+h2)  -  g  ir   =  fv2  -r      -^ 


+  A|>  (2) 


Equations  (1)  and  (2)  are  the  radial  momentum  equations 
in  the  upper  and  lower  layers  respectively.   These  include 
the  advection,  pressure  gradient,  coriolis,  stress,  and 
lateral  diffusion  terms. 

S    I 

3vi  .    3vi      -      t  y-x  y  ,  A32Vi  ^o\ 

-^T1-    +    U^ — L  =  -fUi  +  r +  A-r — ^  (3) 

9t      ldx  *      phi       9x2 


I    B 

|Vi  +  U2|Y^  =  _fu2  +  I  y-T  y  +  A|iv^  (4) 

9t       9x  ph2       9xz 
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These  are  the  tangential  momentum  equations,  comprised  of  the 
same  terms  as  the  radial  momentum  equations,  excepting  the 
pressure  gradient  forces,  which  due  to  symmetry  are  zero. 
The  stress  in  the  momentum  equations  is  defined  by 

S 

t    =  -k  t   sin  a  (5a) 

a 

S 

x    =   k  t   cos  a  (5b) 

a 

where  k  is  the  partitioning  constant  which  determines  how 
much  of  the  surface  stress  will  be  used  for  current  genera- 
tion, and  a  is  the  atmospheric  inflow  angle.   Following 
O'Brien  and  Hurlburt  (1972)  the  stress  components  of  the 
interface  are  defined  as 

I 

x  x  =  p  cq(ui-u2 )  (6a) 

w 

I 

t  y  =  Pwcq(v!-v2)  (6b) 

where    q   =  (qi+q2)/2 

i 
and     q.  =  (u?+v?)2 

The  constant  c  is  the  drag  coefficient  at  the  interface,  and 
is  chosen  to  be  3*10" ^  also  following  O'Brien  and  Hurlburt. 

oil  P 

?Ll   +  9T,ut  _  V^    AM         _  VQS  =  (g) 

at     8x       hi     h,us  Lh}         pch  K    J 

1  i  p  i 
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Equations  (7)  and  (8),  which  describe  the  temperature  and 
depth  of  the  mixed  layer,  are  the  same  as  used  by  Fraim 
(1973)  and  Trapnell  (1974).   The  term  M  is  the  contribution 
to  depth  change  due  to  entrainment  mixing.   The  increase  in 
potential  energy  of  the  column  through  entrainment  at  the 
interface  is  equal  to  the  excess  of  mechanical  generation 
(G)  of  kinetic  energy  over  dissipation  (D),  plus  a  contribu- 
tion due  to  convective  overturning  if  the  column  is  losing 
heat  (QF+Qo  >  0).   A  is  called  the  entrainment  parameter 
which  is  1  when  M  is  positive  and  zero  when  M  is  negative. 
Thus  mixing  across  the  base  of  the  mixed  layer  occurs  only 
when  M  is  positive,  and  therefore  acting  to  deepen  and  cool 
the  mixed  layer.   If  M  were  permitted  to  be  negative  it  would 
simulate  heat  extraction  from  the  lower  layer,  to  warm  the 
already  warmer  upper  layer.   As  this  model  depicts  a  wind 
dominated  region,  M  is  always  1. 


|h^  +  |h_2U2.  =  _M  (   9) 

dt      o  X 


If  +  f^  +  4?  -  o  ™ 


Equations  (9)  and  (10)  are  the  equivalent  equations  to 
(3)  and  (4)  but  in  the  lower  layer.   Mixing  does  not  contrib- 
ute directly  to  the  change  in  temperature  at  that  level, 
therefore  only  the  advection  and  vertical  motion  terms  are 
present.   It  can  be  shown  that  with  the  assumed  linear  dis- 
tribution of  currents  and  temperatures  in  the  lower  layer, 
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that  the  temperature  advected  by  the  average  current  compo- 
nent is  that  from  a  level  2/3  of  the  distance  from  the  bot- 
tom of  the  lower  layer. 

w1=  -h2|^  (11) 


Th  =  2T2  -  Tb  (12) 


Equations  (11)  and  (12)  are  diagnostic  equations  used  to 
close  the  equation  set.   Equation  (11)  was  derived  by  inte- 
grating the  continuity  equation  in  the  lower  layer  and  using 
the  boundary  conditions  of  no  mass  flux  through  the  bottom 
of  that  layer.   Equation  (12)  utilizes  the  fixed  bottom  tem- 
perature, the  linear  lapse  rate  assumption,  and  the  newly 
computed  values  of  the  depth  and  average  temperature  of  the 
lower  layer  to  compute  a  new  value  of  temperature  at  the  top 

of  the  thermocline,  T,  . 

'   h 

The  boundary  conditions  allowed  no  heat  or  momentum  flux 
through  the  ends  or  bottom  of  the  model.   At  the  surface  a 
momentum  flux  from  the  atmosphere  to  the  ocean  was  applied. 
A  heat  flux  from  the  ocean  to  the  atmosphere  also  occurred 
at  the  surface.   The  model  was  considered  to  be  a  free  sur- 
face model  because  the  height  of  the  ocean  surface  varied 
with  respect  to  the  fixed  bottom  depth. 

As  in  O'Brien  and  Hurlburt  the  solutions  to  the  equa- 
tions contain  external  and  internal  gravity  waves.   There- 
fore the  pressure  gradient  terms  in  the  radial  momentum 
equations  were  treated  implicitly  to  allow  a  longer  time 
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step  to  be  taken  and  still  retain  these  waves.   Treatment  of 
these  equations  is  shown  in  Appendix  A.   The  remainder  of 
the  equations  were  evaluated  using  centered  differencing  in 
time  and  space.   Stability  problems  will  occur  due  to  using 
centered  differences  in  both  time  and  space  in  the  diffusion 
term,  therefore  time  lag  was  used  in  this  term. 

The  atmospheric  parameters  were  computed  at  100  grid 
points  with  a  three  kilometer  spacing.   This  atmospheric 
forcing  was  moved  at  a  constant  speed  across  this  grid  in 
order  to  simulate  the  motion  of  a  storm.   The  atmospheric 
forcing  was  applied  to  the  ocean  model  as  shown  in  figure  1. 
Between  and  including  C  and  D  the  atmospheric  forcing  param- 
eters were  applied  directly,  from  D  to  E  the  values  at  D 
were  applied.   The  values  from  C  to  D  in  the  atmosphere  were 
applied  to  the  points  from  C  to  B  in  the  ocean,  although  the 
sign  was  changed  on  the  stress  terms  for  symmetry.   From 
point  B  to  point  A  the  values  at  point  B  were  applied. 
Points  B,  C,  and  D  move  from  right  to  left  to  simulate  the 
motion  of  the  storm. 

As  in  Trapnell's  model  the  surface  stress  and  heating 
functions  were  not  initially  applied  at  full  value.   To  per- 
mit the  initially  quiet  ocean  to  adjust,  the  forcing  was 
increased  exponentially,  requiring  about  nine  hours  for  ap- 
plication of  95%  of  the  full  force  of  the  storm.   It  was 
found  that  both  the  atmospheric  and  oceanic  models  remained 
stable  for  time  steps  of  900  seconds. 
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B.   ATMOSPHERIC  MODEL 

The  atmospheric  hurricane  model  used  to  drive  the  ocean 
in  these  experiments  is  the  two-layer,  time-dependent  model 
of  Elsberry,  Pearson,  and  Corgnati  (1974)  which  was  also 
used  in  the  moving  model  by  Trapnell  (1974).   The  hurricane 
was  considered  to  be  mature  so  that  an  axially-symmetric 
storm  could  be  assumed.   The  maximum  winds  were  permitted  to 
intensify  to  50  m  sec-1  and  then  held  fixed.   This  allowed 
the  changes  in  the  ocean  structure  to  be  isolated  from 
changes  in  the  atmospheric  model. 

The  atmospheric  model  was  integrated  simultaneously  dur- 
ing these  experiments  using  the  new  sea-surface  temperatures 
computed  by  the  ocean  model.   New  surface  stresses  and  sur- 
face heating  values  were  also  applied  to  the  ocean.   This 
constituted  the  interaction  between  the  ocean  and  the  atmos- 
phere. 

This  model  consisted  of  a  surface  inflow  layer  and  an 
upper  outflow  layer.   The  tangential  wind  flow  at  the  surface 
outside  the  radius  of  maximum  winds  was  given  by 

i 
Vfir2  =  const 

and  inside  the  radius  of  maximum  winds  by 

VQr_1  =  const. 

These  equations  were  empirically  derived  by  Riehl  (1963)  and 
supported  by  observations  by  Gray  and  Shea  (1973).   The 
radius  of  maximum  wind  was  assumed  to  coincide  with  the 
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eyewall  radius.   Normally  in  this  model  the  eyewall  radius 
is  defined  to  be  that  radius  at  which  the  average  of  five 
adjacent  values  of  equivalent  potential  temperature  was  a 
maximum.   However,  to  isolate  the  changes  in  the  ocean,  the 
eyewall  radius  was  fixed  at  30  km  after  a  short  period  of 
adjustment.   The  boundary  layer  in  the  atmospheric  model  was 
a  two-layer,  baroclinic  model  developed  by  Cardone  (1969). 
The  heat  and  momentum  fluxes  at  the  interface  were  computed 
in  this  model.   These  were  a  function  of  the  stability  of 
the  boundary  layer  and  the  wind  at  the  top  of  this  layer. 
To  obtain  more  information  regarding  the  atmospheric  model 
the  paper  by  Elsberry,  Pearson,  and  Corgnati  (1974)  should 
be  consulted. 
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III.   RESULTS 

A.   EFFECTS  OF  TRANSLATION  SPEED 

The  simulated  storm  was  moved  at  various  speeds,  and 
compared  with  the  stationary  case  to  observe  the  changes  in 
the  behavior  of  the  ocean  due  to  storm  translation  speed. 
Many  of  the  model  results  are  consistent  with  previous 
theories  and  observations,  as  discussed  earlier.   The  be- 
havior of  selected  ocean  parameters  is  displayed  in  figures 
3  through  6  and  10  through  13.   The  parameters  displayed  in 
these  figures  are: 

a)  the  temperature  at  the  top  of  the  thermocline; 

b)  the  temperature  of  the  mixed  layer,  and  therefore 
the  ocean  surface: 

c)  the  depth  of  the  mixed  layer; 

d)  the  vertical  velocity  as  derived  from  the  continuity 
equation  in  the  lower  layer; 

e)  the  upper-level  current  component  directed  along  the 
path  of  the  storm  (positive  values  indicate  a  cur- 
rent toward  the  right;  also  note  that  the  baseline 
is  at  - . 1  m  sec- l ) ; 

f)  the  tangential  wind  velocity  at  the  top  of  the  atmos- 
pheric boundary  layer  (the  stress  outside  the  storm 
is  held  at  the  value  present  at  the  outer  edge  of 
the  storm,  rather  than  being  zero  as  implied  by 
this  graph). 
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In  the  standard  model  half  of  the  surface  stress  was 
used  to  produce  currents,  and  half  was  used  for  energy  input 
to  entrainment  mixing.   The  total  surface  stress  had  a  maxi- 
mum value  of  2.6  nt  m~2  (26  dynes  cm-2).   The  mixed  layer 
had  an  initial  depth  of  30  m  and  temperature  of  30°C,  while 
the  lower  layer,  the  thermocline,  extended  100  m  below  the 
mixed  layer  and  terminated  with  a  fixed  bottom  temperature 
of  20°C.   This  initial  state  is  depicted  in  figure  2.   As 
can  be  seen  from  figure  2,  a  small  ledge  in  the  temperature 
appears  at  the  bottom  of  the  mixed  layer.   From  Eq.  (7)  it 
is  seen  that  if  T   is  equal  to  T,  the  term  M  becomes  infi- 
nite.  To  maintain  a  slightly  stable  layer  at  the  interface, 
this  difference  is  initialized  at  .1  degree  and  never  al- 
lowed to  become  smaller  in  magnitude.   However,  as  this  tem- 
perature difference  immediately  increases  in  response  to  the 
wind  mixing,  the  condition  was  not  forced  following  initial- 
ization . 

1.   Stationary  Storm  Model 

The  results  after  integration  for  36  hours  with  the 
translation  speed  equal  to  zero  are  displayed  in  figure  3. 
The  wind  profile  is  shown  in  figure  3f.   Note  that  this  does 
not  depict  the  wind  outside  the  storm  area.   The  vertical 
motion  at  the  bottom  of  the  mixed  layer  is  displayed  in  fig- 
ure 3d.   Combining  this  with  the  radial  currents  in  the 
upper  layer,  as  shown  in  figure  3e,  we  see  there  are  two 
circulation  cells  on  each  side  of  the  storm  center.   The 
first  cell  has  its  maximum  upwelling  at  the  inside  edge  of 
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L    =    30  m 


T      =   T.    =    30.00°C 
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Tu    =    29.90°C 
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Figure  2.   Initial  Ocean  Structure 
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Figure   3.       36   Hour  Predicted  Values   for   the   Stationary 
Storm  Case. 
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Figure  3.   36  Hour  Predicted  Values  for  the  Stationary 
(Cont.)    Storm  Case. 
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the  radius  of  maximum  winds,  RMW,  and  downwelling  at  the 
outer  edge  of  the  RMW.   A  second  cell  of  about  half  the 
magnitude  is  located  at  about  90  km  from  the  center.   The 
maximum  depth  of  the  mixed  layer  is  observed  at  the  outer 
edge  of  the  RMW,  rather  than  coincident  with  it  as  might  be 
expected.   At  the  RMW  the  deepening  due  to  mixing  is  offset 
by  upwelling. 

Referring  to  figures  3b  and  3d  the  minimum  tempera- 
ture of  the  mixed  layer  does  not  correspond  to  the  maximum 
upwelling  region,  or  from  figure  3c,  to  the  maximum  depth  of 
the  mixed  layer.   Instead  the  lowest  temperatures  are  di- 
rectly under  the  RMW,  indicating  that  both  upwelling  and 
mixing  contribute  to  cooling  of  the  upper  layer.   The  tem- 
perature at  the  top  of  the  thermocline  is  shown  in  figure 
3a.   Deepening  of  the  mixed  layer,  and  the  resultant  reduc- 
tion in  thickness  of  the  lower  layer,  would  be  the  greatest 
factor  in  the  lower  temperature  at  the  top  of  the  lower  lay- 
er.  The  reduced  thickness  would  lower  the  temperature  di- 
rectly due  to  the  lower  position  on  the  thermal  lapse  rate. 
The  water  being  advected  into  this  region  would  also  be 
colder,  thus  contributing  to  lower  temperatures  (see  the 
discussion  of  T2  in  the  model  equation  (10)). 

Table  I  compares  Trapnell's  Ekman  current  model  with 
this  dynamic  model.   The  results  are  for  a  stationary  storm 
case  after  being  integrated  18  hours.   The  major  difference 
is  in  the  minimum  mixed  layer  depth.   This  could  be  due  to 
the  slower  reaction  of  the  dynamic  model,  as  this  value  was 
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about  24  m  after  27  hours.   Another  explanation  could  be  the 
different  method  of  partitioning  the  surface  stress.   From 
this  comparison  the  results  from  a  model  using  currents  in 
Ekman  balance  appear  to  compare  favorably  with  the  dynamic 
model. 

2.   Moving  Storm  Model 

The  results  of  moving  the  storm  at  3.2  kt  for  72 
hours  are  shown  in  figure  4.   Comparing  figure  4b  with  fig- 
ure 3b,  it  is  seen  that  a  minimum  mixed-layer  temperature  of 
28.3°C  was  observed  as  compared  to  27.3°C  associated  with 
the  stationary  storm.   Evidence  of  waves  on  the  interface  at 
the  bottom  of  the  mixed  layer  is  shown  in  figure  4c.   The 
wavelength  is  about  180  km.   These  waves  are  more  noticeable 
in  the  vertical  and  upper  layer  radial  velocities  in  figures 
4d  and  4e. 

The  comparison  of  Trapnell's  Ekman  model  moving  at 
3.39  kt  and  the  dynamic  model  moving  at  3.2  kt  is  shown  in 
Table  II.   The  main  difference  is  the  momentum  terms,  how- 
ever overall  these  values  of  maximums  and  minimums  do  not 
differ  greatly.   The  problem  with  the  Ekman  model  is  that  it 
does  not  contain  the  dynamics  necessary  to  produce  the  hori- 
zontal distribution  of  variables  that  are  present  in  figure 
4.   It  appears  that  for  sea-surface  temperature  prediction 
the  Ekman  balance  model  would  produce  comparable  results  to 
the  more  complicated  dynamic  model.   However  for  a  predic- 
tion of  the  internal  structure  of  the  ocean,  (e.g.  mixed- 
layer  depth,  slope  of  the  thermocline,  etc.)  the  dynamic 
model  is  necessa.ry. 
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Figure  4.   72  Hour  Predicted  Values  With  the  Storm 
Moving  at  3.2  kt . 
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The  translation  speed  of  the  storm  was  increased 
from  3.2  kt  to  6.5  kt .   In  figure  5  the  results  of  inte- 
grating this  model  36  hours  is  presented.   This  can  be  com- 
pared with  figure  4  as  the  same  distance  has  been  traveled 
by  the  storm  but  in  one-half  the  time.   Note  the  longer 
wavelength  in  the  faster  storm  and  also  the  lag  of  the  up- 
welling  region  behind  the  storm  center  is  much  longer  in 
this  case.   Notice  also,  comparing  4b  and  5b,  that  less 
cooling  has  occurred  due  to  the  passage  of  the  faster  storm 
with  a  minimum  mixed-layer  temperature  of  28.9°C  compared 
with  28.3°C  in  the  slower  storm.   At  the  top  of  the  thermo- 
cline  the  same  results  occur  with  26.5°C  for  the  faster 
storm  and  25.0°C  for  the  slower  storm.   This  case  was  also 
integrated  to  72  hours  and  these  results  are  presented  in 
figure  6.   This  will  be  considered  the  standard  case  for  the 
following  discussion.   Notice  that  there  is  little  differ- 
ence compared  with  the  36  hour  case  in  figure  5,  except  for 
the  longer  wake  behind  the  storm.   Temperatures  were  only 
slightly  lower,  by  .3°C  for  the  top  of  the  thermocline  and 
.1°C  for  the  mixed  layer.   The  wavelength  and  magnitude  of 
the  variations  in  mixed-layer  depth,  vertical  velocity, 
radial  velocity  and  temperature  at  the  top  of  the  thermo- 
cline remain  essentially  the  same.   This  wavelength  is  ap- 
proximately 360  km  as  compared  to  180  km  in  the  slower  storm, 

These  results  may  be  compared  with  the  data  from 
Sheets  (1974).   In  the  observations  the  temperature  at  the 
interface  is  higher  for  greater  mixed-layer  depths,  MLD ,  and 
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Figure  5.   36  Hour  Predicted  Values  With  the  Storm 
Moving  at  6.5  kt . 
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Figure  6.   72  Hour  Predicted  Values  With  the  Storm 
Moving  at  6.5  kt. 


39 


O 

d) 

m 

_J 

CL 

. — . 

o 

CL 

1 — 1 

~1 

\-- 

n 

lu 

cr 

m 

rx_ 

ii  i 

\ 

* 

z? 

ZZ 

U 

T 

CL 

. 

i-u 

_J 

Li 

LiJ 

LO 

CO 

L> 

- 

e) 


CL 


CO 


1 — 1 

L4_. 

r\l 

□ 

Ul 

CL 

\ 

r? 

CL 

ZZ 

z: 

CL 

* 

LiJ 

.J 

CJ 

LlJ 

^r 

CO 

Z> 

1 

f) 

CJ 


CL 

CO 

i — i 

\ 

h- 

z: 

z: 

>— ' 

t  i.i 

LT 

C?i 

z: 

z: 

CL 

\—  i 

f— 

L? 

ID 
<X1 


O 


G 


309 


BOO     900 
DISTRNCE  I  KM- 


1200 


1500 


Figure  6.   72  Hour  Predicted  Values  With  the  Storm 
(Cont.)    Moving  at  6.5  kt . 


40 


lower  for  shallow  MLD's,  which  implies  that  the  waves  are 
due  to  upwelling  and  downwelling  rather  than  mixing.   In 
figure  7,  derived  from  the  data  of  the  standard  case  we  see 
that  this  same  relationship  occurs.   However,  in  Sheets' 
data  these  temperature  fluctuations  also  appear  within  the 
mixed  layer,  whereas  these  fluctuations  are  not  present  in 
figure  6b.   Also  in  the  observations  the  wavelength  is  on 
the  order  of  35  to  40  km  which  does  not  agree  with  the  model 
wavelength  of  about  360  km.   Geisler's  theory  predicts  the 
longer  wavelengths,  about  185  km  for  a  5.2  knot  storm.   More 
data  with  different  storm  speeds  are  necessary  to  compare 
with  the  theoretical  and  model  predictions  of  wavelengths. 

The  vector  contributions  to  acceleration  of  currents 
in  both  the  upper  and  lower  layers  are  presented  in  figure 
8.   These  values  are  computed  at  the  center  of  the  region. 
The  sum  of  the  vectors  is  proportional  to  the  instantaneous 
acceleration  at  that  point.   The  first  set  of  vectors,  at  27 
hours,  represents  the  ocean  condition  prior  to  passage  of 
the  storm.   At  this  time  the  forces  are  very  nearly  in  bal- 
ance, except  for  the  lower-level  current  component  normal  to 
the  storm  track  (v2).   The  storm  passes  this  point  after  37 
hours.   However,  at  36  hours  the  stresses  have  been  reduced 
to  a  very  small  value  compared  with  the  pressure  gradient 
and  coriolis  terms.   It  is  seen  that  this  sudden  reduction, 
and  subsequent  reversal  of  the  stress  pattern,  produces 
strong  inertial  waves  in  the  currents.   It  should  also  be 
noted  that  the  phase  of  the  pressure  gradient  term  also 
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contributed  to  the  imbalance  of  accelerations.   At  45  hours 
this  term  and  the  coriolis  term  are  in  the  same  direction, 
thus  contributing  to  a  large  imbalance.   This  condition 
would  not  occur  without  the  lag  of  the  pressure  gradient 
term  behind  the  center  of  the  storm.   After  45  hours  it  is 
seen  that  the  coriolis  acceleration  dominates  the  solutions, 
exchanging  kinetic  energy  between  the  current  components 
directed  parallel  (u)  and  normal  (v)  to  the  storm  track.   To 
better  illustrate  this,  a  hodograph  at  point  250  of  the  cur- 
rent direction  and  magnitude  as  a  function  of  time  is  shown 
in  figure  9a.   Before  the  storm  passes,  the  winds  are  from 
the  top  of  the  figure  producing  a  stress  which  creates  cur- 
rents from  that  direction,  and  due  to  coriolis  there  is  a 
current  away  from  the  storm  center.   This  is  the  same  behav- 
ior noted  in  the  stationary  case.   In  the  lower  layer  the 
pressure  gradient  term  produces  currents  toward  the  center 
with  the  coriolis  acceleration  producing  a  current  component 
perpendicular  to  the  path  of  the  storm.   The  storm  passes 
from  the  right  at  about  37  hours.   The  inertial  effects  are 
then  seen  to  begin,  with  the  current  direction  reversed  be- 
tween 36  and  45  hours.   Note  that  the  inertial  period  for 
the  position  of  the  storm,  25°N,  is  28.34  hours,  and  observe 
that  from  t=45  hours  to  t=72  hours,  almost  one  rotation  of 
the  current  direction  has  occurred  in  27  hours.   Therefore 
we  see  that  inertial  effects  are  very  strong,  and  indeed 
continue  to  dominate  after  the  passage  of  the  storm  center. 
Due  to  this  relation  to  the  inertial  period,  and  because  the 
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waves  move  with  the  same  speed  and  direction  as  the  storm, 
it  becomes  apparent  that  the  wavelength  should  be  greater 
for  faster  moving  storms.   This  relationship  is  linear  in 
the  cases  studied  here,  that  is  when  the  storm  movement 
speed  is  doubled  the  wavelength  is  doubled. 

A  hodograph  for  the  lower  layer  at  this  same  point 
and  for  the  same  storm  speed  is  presented  in  figure  9b.   The 
currents  in  the  lower  layer  have  the  same  features  as  the 
upper  layer  except  the  magnitudes  are  smaller  and  are  180° 
out  of  phase  with  the  upper  layer.   Remembering  that  the 
vertical  motion  is  a  function  of  the  first  derivative  of  the 
currents  in  the  lower  layer,  it  is  seen  that  vertical  cur- 
rents are  90°  out  of  phase  with  the  fluctuations  of  lower 
layer  current  parallel  to  the  storm  track.   This  upwelling 
and  downwelling  produces  the  waves  at  the  bottom  of  the 
mixed  layer,  therefore  their  wavelength  is  also  the  same  as 
the  fluctuations  in  currents. 

B.   VARIATION  OF  OCEAN  PARAMETERS 

To  check  the  sensitivity  of  the  model,  the  standard  run 
was  repeated  while  varying  individual  parameters.   In  this 
section  the  effects  of  increased  stability  in  the  thermo- 
cline  layer,  increased  initial  mixed-layer  depth,  increased 
stress  for  currents,  and  the  elimination  of  the  mixing  term 
will  be  discussed. 

Changing  the  temperature  of  the  undisturbed  layer,  and 
hence  the  temperature  at  the  bottom  of  the  thermocline, 
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changes  the  stability  of  that  layer.   In  this  experiment, 
shown  in  figure  10,  the  temperature  was  lowered  to  15°C, 
thus  increasing  the  stability.   The  major  differences  were  a 
slightly  cooler  mixed  layer,  about  .1°C  cooler  than  the 
standard  model,  and  a  significantly  cooler  temperature  at 
the  top  of  the  thermocline,  1.8°C  cooler  than  the  standard 
case.   The  reason  for  the  relatively  small  change  in  the 
mixed-layer  temperature  was  that  it  required  more  energy  to 
lift  the  cooler  subsurface  water  into  the  mixed  layer  through 
mixing.   The  significant  reduction  in  the  temperature  at  the 
top  of  the  thermocline  is  due  to  the  cooler  water  being  ad- 
vected  into  the  region  due  to  the  steeper  lapse  rate  (refer 
to  equation  (10)).   The  currents  and  depths  remained  essen- 
tially the  same. 

The  next  variation  from  the  standard  case  was  to  ini- 
tialize the  mixed-layer  depth  at  50  m.   These  results  are 
shown  in  figure  11.   Due  to  the  larger  amount  of  mass  pres- 
ent within  the  layer  due  to  the  increased  mixed-layer  depth, 
the  velocities  in  the  upper  layer  were  smaller,  with  a  maxi- 
mum of  . 36  m  sec-1  compared  to  . 48  m  sec-1  for  the  standard 
case.   Also,  the  mixed  layer  was  .5°C  warmer  in  this  case. 
Again  the  wave  structure  remained  unchanged  although  the 
overall  depth  of  the  mixed  layer  was  increased. 

The  results  of  changing  the  partitioning  of  surface 
stress  to  70%  for  current  production  and  30%  for  mixing  are 
shown  in  figure  12.   As  expected,  the  currents  were  stronger 
in  this  case  than  in  the  standard  case  of  50%  partitioning 
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Figure  10.   72  Hour  Predicted  Values  With  the  Storm  Moving 
(Cont.)     at  6.5  kt  except  with  the  Initial  Tb  =  15°C. 
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Figure    11.       72   Hour   Predicted   Values   With   the   Storm  Moving 
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shown  in  figure  6.   Comparing  figure  6b  with  12b  it  is  seen 
that  the  surface  temperature  is  lower  when  the  mixing  is 
stronger.   Figures  6a  and  12a  show  the  same  to  be  true  for 
the  temperature  at  the  top  of  the  thermocline.   This  result 
and  the  deeper  mixed  layer  in  6c  compared  to  12c,  show  that 
mixing  is  the  most  important  contributor  to  deepening  of  the 
mixed  layer.   A  close  comparison  of  the  waves  at  the  bottom 
of  the  mixed  layer,  however,  shows  that  the  amplitude  is 
increased  by  about  6  m.   It  therefore  appears  that  the  am- 
plitude of  these  waves  is  a  function  of  the  magnitude  of  the 
currents. 

A  run  was  also  made  with  the  mixing  term  set  equal  to 
zero.   Comparing  figure  13  with  those  of  the  standard  case, 
figure  6,  it  is  seen  that  in  the  case  without  mixing  the 
surface  and  thermocline  temperatures  have  only  been  slightly 
affected.   The  mixed  layer  has  not  deepened  significantly 
although  the  waves  were  present  with  about  the  same  ampli- 
tude.  The  currents  were  also  stronger  even  though  the 
stress  component  partitioned  to  current  production  is  the 
same.   This  effect  is  due  to  the  mixed  layer  remaining 
shallow,  and  therefore  the  magnitude  of  the  depth-averaged 
currents  are  larger. 
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IV.   CONCLUSIONS 

The  results  shown  in  the  previous  section  agree  quite 
well  with  previous  theories  and  observations  regarding  the 
behavior  of  the  ocean  beneath  a  moving  hurricane.   The  cool- 
ing of  the  surface  temperature,  as  observed  by  Leipper 
(1967),  Black  and  Mallinger  (1972),  and  Sheets  (1974),  was 
produced  by  this  model.   Entrainment  mixing  was  the  dominate 
factor  in  cooling  of  the  surface  temperature  and  deepening 
of  the  mixed  layer.   Upwelling  maxima  and  minima  as  pre- 
dicted by  Geisler  (1970)  were  also  produced  in  the  wake  of 
the  moving  storm.   Increased  translation  speed  produced  less 
cooling  behind  the  storm;  this  agreed  with  the  data  of  Black 
and  Mallinger.   In  comparison  with  Trapnell  (1974)  this 
model  appears  to  be  quite  consistent  thermodynamically . 
However,  as  stated  previously  the  earlier  model  did  not  con- 
tain the  dynamics  to  produce  the  horizontal  structure  pro- 
duced by  the  present  model. 

The  primary  discrepancy  was  the  wavelengths  of  the 
waves  at  the  bottom  of  the  mixed  layer.   Sheets  shows  values 
on  the  order  of  50  km  for  a  fast  moving  storm  whereas  this 
model  shows  wavelengths  on  the  order  of  360  km  for  storms 
moving  at  6.5  kt .   This  difference  has  not  been  resolved. 
In  this  model  these  waves  appear  as  a  result  of  inertial 
waves  produced  by  the  sudden  change  in  stress  as  the  storm 
passes  and  due  to  the  lag  of  the  pressure  gradient  force 
behind  the  storm  center. 
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In  addition  to  resolving  the  difference  in  wavelength 
mentioned  above,  other  work  is  necessary  on  this  model.   The 
method  of  partitioning  stress  for  mixing  and  current  produc- 
tion needs  to  be  studied.   Turbulent  momentum  transfer  at 
the  interface  should  be  included.   Sensitivity  tests  should 
be  run  on  the  drag  coefficient  at  the  interface  and  the  dif- 
fusion constant. 

The  most  valuable  refinement  would  be  to  extend  this 
model  to  three  dimensions.   This  would  eliminate  the  most 
restrictive  assumption  in  this  model,  axial  symmetry.   To 
make  this  modification,  and  maintain  the  present  resolution 
in  the  most  active  region,  near  the  hurricane  center,  would 
require  some  type  of  nested  grid  arrangement.   The  nested 
grid  approach  would  be  necessary  due  to  the  present  limita- 
tion of  computers  in  speed  and  size. 
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APPENDIX  A 

SEMI-IMPLICIT  SCHEME  FOR 

MOMENTUM  AND  DEPTH  EQUATIONS 

The  equations  in  finite  difference  form  which  were  used 
in  the  computer  model  are  developed  here.   Using  the  equa- 
tions presented  in  the  discussion  of  the  model  we  proceed  to 
implicit  solutions,  for  u1?  u2,  hlt    and  h2.   Following  the 
procedure  in  O'Brien  and  Hurlburt  (1972) 

S      I 
+A0J-  e  U,  (A-l) 


If  +  H't^  "  -h.'lf-  -*ifr   +  M  ^  <A-2' 


I 

9u2  ,       %     ,u    j-u    \  '  3hj_      9u?    a  xX 

FT  +  ga^(hi+h2 }  -  g  3^  =  -u*^r  +  fv2  +  ^ 

+4^- E  U2  (A~3) 


|*k.  +  H  |Ui  =  _h'iu2.  _    8A.  .  M  =  D  (A_4) 

8t      ^3x       23x      23x  2  v    ' 

The  tendencies  are  now  evaluated  for  a  2At  time  step  and  the 
divergence  terms  are  averaged  between  time  levels  (n-l)At 
and  (n+l)At.   The  subscript  j  is  the  space  index  while  the 
superscript  n  is  the  time  index. 
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u  .  +  [Atg-~—  (hi+h2  )]  •    =  2AtU,  ,+u,  .  -[Atg-~— (hi+h2  )]  • 
ij      b3x        3  13       *J       dx  3 

=    Li  (A-5) 


v,n+lxrAfu   3ui,n+l        OA,nn    ,,n-l    r  A,  „   8ui,n-l    _    T  ,.    c. 

h,  .    +[AtH!7^ — -]  .         =    2AtD.  .+h,  .    -[AtHi-r — -]  .         =    L2  (A-6) 


n+1,  r ,,     d    , ,     , ,     x    .,     )9hi-,n+l        n,,..!!    ,    n-1 
U2j    +[Atg^(h1+h2)-Atg  j^}.         =    2AtU2j+u2j 

-[Atg|^(h1+h2)-Atg,||i]^"1    =    L3  (A-7) 


i_n+l,r..TT  9u2,n+l   ^.A^T^n  ,,n— 1  taj-tt  9u2ln-l  _  T       •  »  <->\ 
h  .  +[AtH2^— -]  .    =  2AtDo  .+h,  .  -[AtH2-r— -]  .    =  L^      (A-8) 


With  homogeneous  boundary  conditions  on  Ui  and  u2 ,  hi  and  h2 
can  be  eliminated.   Substituting  (A-6)  and  (A-8)  into  (A-5) 
and  (A-7)  respectively,  and  dropping  the  n+1  superscript 
the  following  is  obtained. 


u1.+{Atg|^[L2-(AtH1|^)+U-(AtH2|^)]}   =  L,  (A-9) 


u2j  +  {Atg|^[L2-(AtH1|^-)+LI+-(AtH2|^)] 

-Atg'|^[L2-(AtH1|^)]}J  =  L3  (A-10) 


Expanding  (A-9)  gives: 

ull-(At)*g[H1|i(|Hi-)+H1|5(|^)]    -  LI-4tg|i(LI+LO  (A-ll) 


61 


Expanding  (A-10)  gives: 

u  .-(At)2g[H1|-(|^)+H2|-(|^-)]+g,(At)2H1|-(|^-) 

2j     '      3x  3x      dx    dx  9x  3x 

=  L3-Atg[|^(L2+LO]+Atg'|^  (A-12) 

Subtracting  (A-ll)  from  (A-12)  gives: 

u.j-u^VcAt^H.f^fii-)  =  L,-L1+Atg'fii  =e  (A-13) 

Substituting  u  .  from  (A-13)  into  (A-ll)  gives: 

u2r(At)2[(g-g,)H4(|^)]-g(At)2H4(|^) 

=  L3-Atg|^(L2+LO+Atg'|^-  E  d  (A-14) 

For  the  remainder  of  this  derivation  only, the  following 
definitions  are  made. 

a  =    gHi(At)2 

b  =    gH2(At)2 

c  e  g'Hi(At)2 

f  E  l/(Ax)2 

From  (A-13)  we  get 

u  ,-u,  .+0!-^-  =  e  (A-15) 

2J   ij   3x2 

From  (A-14)  we  get 


U2j-<a-C>i^  "  H^  -  d  (A-16) 
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■}2  u  .     --2u  .+U  .     - 

Using  |^t  =      J    (Ax)2 'LJl     E(i-    (A"15)   becomes 


u„  .-u,  .+cf  (u      x1-2u,  .+u    •    -,  )    =    e  (A-17) 


2J       l3 


(u    ,x1-2u    .+u    .    - )   =   e 


1  u2-i~e 

-u  ._,.., +  (2+  -J7)u  .-u,  .  1  =  — ^ —  (A-18) 

ij+1  v    cfy  ij   ij-l    cf 


This  is  the  form  to  apply  the  method  presented  in  Richtmyer 
(1967).   Using  the  symbol  conventions  in  that  text, 

A  =  1 

B  =  (2  +  -%) 
v     cf 


C  =  1 

J   J 


D  =  (u2,-e.)/cf 


As  c  and  f  are  both  positive  the  criterion,  explained  in 
Richtmyer,  of  B>A+C  is  satisfied. 

Having  used  a  guess  for  u  .  above,  a  new  u2 .  using  the 
above  ux  .  can  now  be  computed  from  Eq.  (A-10). 

%J 

-u       ,.,+(2+  rAr)u     .-u     .    ,    =  ■—   +  ^r^(u     .^-,-2u     .+u     .    -)       (A-19) 
2j+l    v         fby    2j       2j-l        fb  b    v    ij+1         ij       ij-l 

As   used    for   Uj 
A   =    1 
B   "    2    +7b 


C   =    1 


d         a-c 
fb  b    ^ij+1    ""i  j    "i  j- 


D  :  +  "-b-^ii+l-^ii^ij-l) 
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Solve  Eq.  (A-18)  and  (A-19)  iteratively  to  obtain  values  for 
Uj  and  U2 . 

The  scheme  in  Richtmyer  is 

u.  =  E.u.^+F.  (A-20) 

J     J  J+l   J 

where 

A.  1 

Ej  =  B.-C.E.  ,  =  B.-E.  "  (A-21) 


and 


D.+F 

Fj    =  B^^T  (A"22) 

3      J-l 


The  boundary  conditions,  u  ,  u2,  EQ,  and  F  ,  are  zero. 

Values  for  hx  .  are  obtained  by  rearranging  Eq .  (A-6) 


i_n+l    t   f.j.11  3ui,n+l  , .    OON 

h  .   =  L2-[AtHi^ — L]  .  (A-23) 


Likewise  Eq.  (A-8)  provides  a  relationship  for  h  . . 

,n+l    T   rA-nT  3u2,n+l  . .    _.. 

h,  .   =  L^-  [AtH2^r— *■]  .  (A-24) 

This  completes  the  semi-implicit  scheme  used  in  computing 
new  values  for  U],  u2 ,  hi  and  h2. 

The  remainder  of  the  equations  are  handled  explicitly, 
using  centered  differences  in  time  and  space. 
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